// DBSCAN constructs a DBSCAN cluster analyzer.
// This implements the DBSCAN algorithm published in Wikipedia's DBSCAN article
// as of October 19, 2010. I'm including the algorithm as published, with my notes,
// at the end of this file. I've used the same variable names in the Javascript code
// to help correlate the code and algorithm/article.
//
// Call the DBSCAN constructor with four or five arguments:
//      an array of data points     D
//      a distance function         dist
//      an epsilon value            eps
//      a minimum cluster size      MinPts
//      an initialization function  init   (optional)
//
// You may update these as you wish before running or rerunning an analysis.
// Run your DBSCAN by calling its run() method.
// Run() accepts optional epsilon and cluster size arguments. 
// You may also set those and other properties directly before running a DBSCAN.
//
// If you supply an init function, the constructor will call it as a last step
// in constructing the DBSCAN analyzer. The run() method will also call it before
// each DBSCAN run.
//
// D may be an array of primitive values or any kind of object. Only dist()
// needs to understand the structure of D elements to compute the distance
// between any two of them, given their indices. Other DBSCAN routines merely 
// treat data points as index values.

function DBSCAN (D, dist, eps, MinPts, init) {
    // Construct a DBSCAN cluster analyzer.
    this.D = D;           // array of data points
    this.dist = dist;     // distance function(i1, i2) given indices of two data points
    if (init)
      this.init = init;   // optional initialization function, to run upon construction
                          //   and before each analysis run
    this.eps = eps;       // neighborhood radius
    this.MinPts = MinPts; // minimum number of points to form cluster
    this.assigned = [];   // cluster assignment 
                          //          >= 0 --> cluster index
                          //            -1 --> noise
                          //     undefined --> we have not yet visited the point
    this.cluster = [];    // array of clusters, each an array of point indices
    // Note that we store cluster assignments redundantly in this.assigned and this.cluster
    // to quickly determine which points are in a cluster and which cluster a point is in.
    // Always update both arrays!
    this.run = dbscanRun;  // run the analysis, optionally with new eps, MinPts values
    this.getNeighbors = dbscanGetNeighbors;   // private
    this.expandCluster = dbscanExpandCluster; // private
    if (this.init) this.init();
}

function dbscanRun(eps, MinPts) {
    if (eps) this.eps = eps;
    if (MinPts) this.MinPts = MinPts;
    this.assigned = new Array(this.D.length);
    this.cluster = new Array();
    if (this.init) this.init();
    for (var P in this.D) {
        if (this.assigned[P] === undefined) {  // if we haven't visited this point already
            console.log('visiting ' + P);
            var N = this.getNeighbors(P);
            // console.log(N.length + ' neighbors: ' + N);
            if (N.length + 1 < this.MinPts) {
                this.assigned[P] = -1;  // noise
                // console.log('noise');
            }
            else {
                var C = this.cluster.length;  // next cluster index
                this.cluster[C] = [];         // new cluster
                // console.log('cluster ' + C);
                this.expandCluster(P, N, C);
            }
        }
    }
    // for (var c in this.cluster) {
    //     console.log('cluster ' + c + ' size ' + this.cluster[c].length + ': ' + this.cluster[c]);
    // }
}

function dbscanGetNeighbors(P) {
    var neighbors = [];
    for (var i in this.D) {
        if (i == P) continue;
        if (this.dist(P, i) <= this.eps)
            neighbors.push(i);
    }
    return neighbors;
}

function dbscanExpandCluster(P, N, C) {
    this.cluster[C].push(P);  // add point index to cluster
    this.assigned[P] = C;     // store cluster index in point assignment
    for (var PP = 0; PP < N.length; PP++) {  // PP means P' -- note P' is indexing N, not D
        // console.log('> ' + N[PP]);
        if (this.assigned[N[PP]] === undefined) {  // P' not yet visited?
            // console.log('> visiting ' + N[PP]);
            var NP = this.getNeighbors(N[PP]);     // NP means N'
            // console.log('> ' + NP.length + ' neighbors: ' + NP);
            if (NP.length + 1 >= this.MinPts) {
                N = N.concat(NP.filter(function(p) { return p != P && N.indexOf(p) == -1 } ));
                // console.log('expanded neighborhood: ' + N);
            }
        }
        if (!(this.assigned[N[PP]] > -1)) {  // P' not yet assigned to a cluster?
            this.cluster[C].push(N[PP]);  // add point index to cluster
            this.assigned[N[PP]] = C;     // store cluster index in point assignment
        }
    }
}

/* DBSCAN algorithm

DBSCAN(D, eps, MinPts)
   C = 0
   for each unvisited point P in dataset D
      mark P as visited
      N = getNeighbors (P, eps)
      if sizeof(N) < MinPts
         mark P as NOISE
      else
         C = next cluster
         expandCluster(P, N, C, eps, MinPts)

expandCluster(P, N, C, eps, MinPts)
   add P to cluster C
   for each point P' in N 
      if P' is not visited
         mark P' as visited
         N' = getNeighbors(P', eps)
         if sizeof(N') >= MinPts
            N = N joined with N'
      if P' is not yet member of any cluster
         add P' to cluster C

= Discussion =

== definitions for directly density-reachable, density-reachable, and density-connected ==

q is directly density-reachable from p if
  q is in p's epsilon neighborhood, and
  p's epsilon neighborhood contains enough points to constitute a cluster

q is density-reachable from p if
  there is a sequence p(1),...,p(n) where
    p(1) = p
    p(n) = q
    each p(i+1) is directly density-reachable from p(i)

The density-reachable relationship is not commutative--q may lie on the fringe
  of a cluster and may not itself have enough neighors to originate a cluster.
Thus we say p and q are density-connected if there is a point o such that 
  p and q are each density-reachable from o.
  p  <-- o --> q


== cluster definition in terms of density-connected points ==

A cluster is a subset of the points in a database where
  1. All points in the cluster are density-connected.
  2. Any point density-connected to a point in the cluster is part of the cluster.


== algorithm restated in plain language ==

DBSCAN parameters are
  epsilon - the radius, in terms of the distance function dist(), 
            of a neighborhood around each point
  minimum points - minimum number of point in an epsilon neighborhood 
                   to originate or extend a cluster

Select an arbitrary starting point P that we have not yet visited.
Retrieve the points in the epsilon neighborhood around P.
If we have at least the minimum number of points,
  originate a new cluster containing the points from P's epsilon neighborhood.
Otherwise,
  Label the point P as noise (not yet a member of a cluster).
  Later we may find the point in a sufficiently populated epsilon neighborhood of another point and and assign it to that new cluster.

For each point in a new cluster,
  Retrieve the points in its epsilon neighborhood.
  If the number of points is at least the minimum,
    add those points that are not already in the cluster to the cluster.
  Continue visiting points and expanding the cluster until you have visited
    every point in the cluster.

Select a new, unvisited point P and repeat from the top, until we have
  visited every point in the database.

*/
